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Abstract: 

Recent advances in interfacing computational software with experimental hardware have 
expanded the area of application of virtual reality (VR) beyond that of visualization to simulation- 
based training and education. The VR can now be directly interfaced with hardware and physics-based 
computational models to form a sophisticated hybrid simulation system. The integrated simulation 
system developed in litis research demonstrates the motion of a heaving buoy in irregular waves in 
order to demonstrate lite process of wave energy conversion. The VR system consists of three main 
function blocks: A physics-based simulator in MATLAB® , a game studio in Unity 3D, and a hardware 
manipulator in LabVIEW®. These three blocks operate in a decentralized mode and communicate with 
each other as interconnected nodes of an integrated network. Tin's article details the development of 
each of these function blocks and addresses issues associated witlt data transmission and lockstep 
synchronization. The effectiveness of the developed hybrid VR system is demonstrated to conclude the 
presentation. 
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I. INTRODUCTION 

This research follows the Mikropoulos and Natsis [1] notion of virtual reality (VR), who define it as "a 
mosaic of technologies that support the creation of synthetic, highly interactive three-dimensional spatial 
environments that represent real or non-real situations." The dominant technology in VR is visual 
representation, namely virtual environment (VE). VE and VR are quite often interchangeable in the literature. 
VE along with technologies such as immersion, multisensory interaction and haptic devices has enabled VR to 
offer a first-order experience in emulating the natural environment and visualizing abstract ideas. As a result, 
VR has been proven to be effective in increasing knowledge construction and, thus, enhancing the learning 
outcome. 

Many educational institutions have recently incorporated VR in classroom teaching; particularly in the engineering 
and science disciplines such as factory production process [2-3], chemistry [4], electrical [5] and mechanical engineering 
laboratories [6]. Recent advances in software-hardware interfacing have increased the level of the VE sophistication and 
broadened the range of applications of VR to allow the VE to interact with physical models in real time. For instance, Hu 
and Zeigler [7] developed an interface tread to allow for direct interaction between a robot and the VE. Richardson, et al [8] 
incorporated the machine code in the virtual environment to simulate the microcontrollers interactively. Morris et al [9] 
developed a physics-based simulation to model sound and force dynamics in a virtual environment to train physicians for 
temporal bone surgery. Fu et al [10] introduced risk models in the virtual environment for safe navigation training in the 
three gorges reservoir area. In addition to those project- or course-oriented VEs, many educators and IT experts have 
proposed proper frameworks and simulation tools to develop a distributed VE shared by multiple schools [11-13]. 

In order to investigate the effectiveness of using VR as an educational tool, Mikropoulos and Natsis 
recently reviewed 53 research articles on educational virtual environments and summarized their findings in [1] 
in terms of educational context, characteristics and features of each VR, and the associated learning theory. 
Abulrub et al [14] stressed creative learning in an interactive and immersive virtual environment in which the 
students can create and investigate design alternatives. Finally, Lindgren [15] used an expert's view in a virtual 
factory for safety training to prove that changing a person's perspective in a virtual environment benefits 
learning. 
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The goal of this research is to develop a hardware in-the-loop VR simulation system that can be used to 
demonstrate the ocean wave energy harnessing process and the associated technologies. Ocean waves offer an 
attractive green and renewable energy source and have generated considerable interest in research, development 
and test in recent years. For the purpose of this research, the wave energy harnessing is performed by a heaving 
buoy, which is connected to a power take-off device and a generator in order to convert the wave energy first 
into mechanical energy and finally into electricity. The integrated hybrid VR system developed to demonstrate 
the wave energy conversion (WEC) process comprises three major components: A physics-based ocean wave 
simulation and wave-buoy interaction software, a hardware mechanism to convert the heaving buoy motion into 
electricity, and a virtual environment to display the interaction between the buoy and the waves. The WEC 
simulator sends the buoy motion signal simultaneously to both the hardware stand and the VR so that the buoy 
moves in synchrony in real and virtual environments. The challenge faced in this task corresponds to achieving 
a lock-stepped synchronization among the involved three platforms: physics-based simulation, virtual reality, 
and hardware mechanism. 

Drafted in the early 2000s, Unity 3D has been used since as a cross-platform to create many new 3D 
games and applications for mobile devices, desktop PCs, and the web. It is also suitable for creating menu 
interfaces using built-in animation and coding scripts. In this project, Unity 3D is utilized to create the VE based 
on the data provided by the physics-based simulator. In Unity 3D, a game scene is constructed to demonstrate 
the working principle of the WEC process, i.e., absorb the wave energy through the heaving motion of the buoy 
and convert it into mechanical energy. C# scripts are used primarily to control every action in this game scene 
from the sunlight reflection on the surface of the ocean and the propagation of the waves to the motion of the 
floating buoy. At the beginning of each simulation, all the scripts used in the game scene are compiled and later 
utilized to regulate the game object according to the data input from the physics-based simulator. 

The paper is organized as follows: Section 2 presents the wave spectrum model used in the Unity 3D to 
generate ocean waves in Unity 3D. The same wave spectrum is utilized to excite the floating buoy in the 
physics-based simulator. The equation of motion of the heaving buoy that results from the wave-buoy 
interaction and the equation of the corresponding energy absorption are derived in Section 3. The data flow 
among the VR system and the WEC hardware mechanism is described in Section 4. A few examples are given 
in Section 5 to demonstrate the simulation of the wave energy conversion process using the hybrid VR system. 
Conclusions and directions for improvement of the presented system are presented in Section 6. 

H. MODELLING OF OCEAN WAVES 

Considering the stochastic nature of ocean waves, the ocean free surface elevation can be described as 
a random process; for design purposes it is often assumed that it is an ergodic random process, i.e., its statistical 
properties will not change with time and these properties can be deduced from a single, sufficiently long sample 
of the process. Therefore, it can be modeled by a wave spectrum (polychromatic waves), which distributes the 
wave energy over wave frequency and direction. The spectral density can be computed by utilizing existing 
parametric wave spectra, e.g., the JONSWAP (Joint North Sea Wave Project) spectrum [16], which is capable of 
modeling both developing and fully-developed open-sea conditions. In this work, the unified directional 
spectrum, proposed in [17], is utilized in order to model the ocean waves. The spectrum is modeled considering 
the contribution of both short and long wind-driven waves, 

S k (k)=k-\B, + B h ) (1) 

where S k is the wavenumber spectrum, k is the wavenumber, B t is the spectrum curvature for low frequencies, 
and B h is the spectrum curvature for high frequencies. The wavenumber is defined as a function of the wave 
length A: k=2/r/2. Assuming operation in deep water, the relation between wavenumber and wave 
frequency (dispersion relation) is defined as: 

co 2 =hg (2) 

with g being the acceleration of gravity (g = 9.81 m/s 2 ). The spectrum curvature for low frequencies is 
calculated using the following equation: 



www.ijceronline.com Open Access Journal 



A Hybrid Virtual Reality Simulation System for Wave Energy Conversion 



where c is the wave phase speed (c = calk), c p is the wave phase speed at the spectral peak frequency, dp is the 
generalized Phillips -Kitaigorodskii equilibrium range parameter for long waves, and F p is the long-wave side 
effect function. The parameter a p is computed using the following expression [18]: 



where Q is the inverse wave age defined as Q = U 10 /c p . U 10 is the wind speed at a height of 10 m from the 
water surface. The following relationship is proposed in [17] between the fetch x, i.e., the distance over which 
the wind blows with almost constant velocity, Q, and U 10 , 



i2 = 0 - 84 -( tanh (fe) )) 



with X 0 = 2.2T0 4 . In this research, we assume that the WEC system is operating in fully-developed seas, thus, 
with a typical fetch value of 150 km and U 10 that can vary between 2 and 10 m/s, the corresponding values of Q 
are in the range [0.84, 1.11]. The long-wave side effect function is calculated as: 



where k p is the wavenumber at the spectral peak frequency, L PM is the Pierson-Moskowitz shape spectrum [19], 
and J P is the JONSWAP peak enhancement factor [16] given by the following equations: 

L PM =exp{-^(^) 2 } (7) 

y p = Y r , y = 1.7 for 0.84 < Q < 1 and y = 1.7 + 6-log( Q) for 1< Q < 5 

r = exp I - ^ ) [,o- = 0.08 ■ (1 + 4/T 3 ) (8) 



for high frequencies is calculated as: 

B h = 1 -a m jF m (9) 

where a m is the generalized Phillips-Kitaigorodskii equilibrium range parameter for short waves, c m is the 
minimum phase speed at the wavenumber k m associated with the gravity-capillary peak in the curvature 
spectrum, and F m is the short-wave side effect function. Assuming that k m = 370 rad/m [17], then c m = 0.23 
m/s. The short-wave side effect function is calculated as: 



The equilibrium range parameter is calculated assuming a linear relation between a m and the friction 
velocity u* at the water surface: 



where the constant a 0 is set equal to 1.4T0" following the observation in [17] that this value compares 
favorably with the frequency spectrum parameter values reported in [20]. The calculation of u* is performed 
through the log-law equation for boundary layer flows: 
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with z = 10 m, the universal constant k = 0.4, and z 0 is the near-surface aerodynamic roughness length, which is 
estimated from the following semi -empirical equation [21]: 



z 0 =3.7-10- 5 ^(^) 



In order to demonstrate a more realistic ocean environment, short-crested seas, i.e., seas where the wave 
propagation is not unidirectional relative to the wind direction, are also included in the wave simulation. This 
can be accomplished by multiplying eq. (1) by a spreading function to obtain a directional spectrum. The latter 
is the Fourier transform of the auto-covariance of the free-surface elevation, which is a real and even function 
and, thus, the Fourier series expansion of the spreading function must contain only even harmonics. In this 
research, the spreading function proposed in [17] is utilized: 

<P(k, 0)= ^ ■ [1 + A(fc) ■ cos(2<p)] (14) 

where (p is the wave propagation direction relative to the wind direction and A(/c) is the coefficient of the second harmonic 
of the truncated Fourier series expansion of the spreading function <P{k, 0), which is computed as: 

A(fc) = tanh + 4 ■ + 0.13 ■ f- ■ (^r)"] (15) 



Therefore, by combining the omnidirectional spectrum from eq. (1) and the spreading function from eq. 
(14), we obtain the unified directional spectrum [17]: 



[ f^(fc,0)=i-f(fc,0)-5 ft (fc) (16) 

The computation of the wave elevation r){x,y, t) at the location of the buoy can be performed through a 
linear superposition of N discrete wave components, each with wavenumber k n , and M discrete values of the 
angle (p m , i.e., the angle between the wave propagation direction and the wind direction 0 m , 



r\{x,y, t) = ZnLl Zm=i A /V / fcn ,0 m fc n .A/c n A0 m sin(w n t - k n x cos <p m - k n y sin <p m + s nm ) (17) 

where e nm is a random phase angle uniformly distributed in [0, 2k] and the angle 0 m can take values in [-71, k]. 
The computation of the wave elevation can be facilitated by utilizing an inverse Fast Fourier Transform (FFT). 
In order to apply the inverse FFT, eq. (17) is written in the following form: 



■IK 



m k n .Ak n A(p m cos(k n x cos </> m + k n y sin </> m 



k n .Ak n A(p m sin(k n x cos cp m + k n y sin cp m 



In the WEC simulations, we have used 256 components for k n and 0 m , i.e., N = 256 and M = 256 with 
k n in the range [0.04, 4] rad/m. 



III. WAVE-BUOY INTERACTION AND WAVE ENERGY ABSORPTION 

The time-domain simulation of the wave energy conversion process using a hemispherical (half-submerged) buoy 
is described in this section. The incident wave exerts a force on the buoy. This force can be computed by utilizing the model 
for the short-crested, irregular waves, eq. (16), to obtain the spectrum of the wave excitation force, 



S Fw = \H(co)\ 2 -V ki(t> 
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where H(a>) is the force transfer function, which facilitates the computation of the magnitude of the force per 
unit input, i.e., wave force per unit wave amplitude. According to Faltinsen [22], the force transfer function for a 
heaving spherical buoy can be written as: 

H(to) = ^ P9 ^ (20) 

where B h (a)) is the radiation damping coefficient, which for a hemispherical (half-submerged) buoy can be 
obtained from [23], and p is the fluid density. The wave frequency a> is obtained from the corresponding wave 
number from eq. (2). By combining eqs. (18) - (20) and using an inverse FFT, the time series of the incident 
wave excitation force F w is obtained. 

In addition to the force due to the incident wave, the power take-off mechanism that converts the wave 
energy into mechanical energy exerts a force on the buoy, F PT0 . In this work, we assume that this force is 
modeled as a viscous (linear) damper: 

F PT0 (t) = -B PT0 -t(t) (21) 

where the damping coefficient B PT0 is constant and equal to the radiation damping coefficient at the spectral 
peak frequency, B h (oj p ), and <f (t) is the vertical velocity of the buoy. The heaving motion of the buoy results in 
the production of waves radiated away from the buoy, which causes an additional force exerted on the buoy by 
the fluid, the radiation force, F R : 

F R (t) = -m R (oo) ■ ftt) - £jc (t - TXO)dT (22) 

where m^oo) is the added mass of the buoy at infinite frequency, <f(t) is the vertical acceleration of the buoy, 
and K(t) is the radiation force kernel. For a heaving hemisphere, m R (co) is equal to approximately half the mass 
of the buoy, m b [23]. The radiation force kernel is computed as: 

K{t) = m b \ /"flhOO ■ cos(wt) da> (23) 

The convolution integral in eq. (22) is computed numerically using the trapezoidal rule. Finally, the net 
force between buoyancy and weight of the buoy generated due to the displacement £(t) of the buoy from its 
equilibrium position corresponds to the hydrostatic force, F HS : 

F H s(t) = ~pgA w at) (24) 

where A w is the cross-sectional area of the buoy at the free-surface plane. For a floating buoy with a single 
degree of freedom (heave), the rigid body equation of motion is [24]: 

= F w (t) + F PT0 (t) + F R (t) + F HS (t) (25) 

Equation (25) combined with eqs. (21), (22), and (24) result in the following integro -differential equation: 

(m + m R (oo)X(t) + £jc (t - TX(x)dT + B PT0 ■ at) + pgA w Ut) = F w (26) 

This equation is solved using a 4 th order Runge-Kutta solver with a fixed time step of 0.05s. The 
average absorbed power over a time period Tcan be calculated using the following equation: 

Pave =7/ o r ^™(t)-C(t)dt (27) 

If the instantaneous absorbed power is measured at N discrete time intervals over a time period T, then 
eq. (27) can be rewritten in the following form to compute its root mean squared value: 

P m , rms = ^^F PT0 (t n )x(t n ) = ^Z^B PT0 (x(t n )y (28) 
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The physics-based simulator used in the VR system is based on the aforementioned equations for the 
ocean waves and the wave-buoy interaction and it was developed in MATLAB®. The wind speed can be 
configured to create different settings of the ocean. 

IV. SOFTWARE INTEGRATION AND INTERFACE WITH WEC HARDWARE 
MECHANISM 

The integrated VR system was developed around the three main blocks shown in Fig. l.a: The physics- 
based simulator in MATLAB®, the game engine in Unity 3D, and the hardware controller in LabVIEW®. 
Whereas Unity 3D helps to visualize the waves and the buoy in real time, LabVIEW is responsible for 
controlling the hardware (stepper motor). These three blocks work in a decentralized mode and communicate 
with each other as nodes in a network. The entire WEC system is show in Fig. 1 .b. 



MATLAB 


User Datagram Protocol 


LabVIEW 

Hardware Manipulator 


Physics-based Simulator 






Figure 1 .b The Wave Energy Converter 

In order to capture the motion of the buoy from the physics-based simulator, the data in MATLAB® is 
streamed in a lock step synchronously to both Unity 3D and LabVIEW®. Once the time step and other 
simulation parameters have been selected, the simulation code in MATLAB® generates the motion profile of the 
buoy at a single time interval. The User Datagram Protocol (UDP) is then employed to transfer the buoy motion 
profile at every time interval to the other two blocks: The hardware controller to actuate the WEC mechanism 
motor and the game engine to visualize the simulation. UDP is the networking protocol that grants the 
applications direct access to the data transport service. The real time synchronization is achieved by imposing 
the MATLAB® simulator to send buoy's position and velocity simultaneously to Unity 3D and the hardware 
controller at the interval of 0.1 second clock time. The buoy's motion profile is obtained by solving the buoy 
equation, eq. (26) numerically. 
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The simulation time used in the ODE solver is in general shorter than the real clock time. Once the 
simulator finds the buoy's profile at 0.1 seconds, it will pause temporally until the clock time also reaches 0.1 
second interval and then send out the data wirelessly. The Unity 3D spends its simulation time to generate 
random wave patterns and render their image. The duration of the simulation time for wave generation depends 
on the size of the random numbers associated with wavenumber k_n, and the angle (j)_m for wave generation 
and while the rendering time depends on the grid size. The component number is set to be 256 for both 
wavenumber and angle and the grid size is 300x300. The Unity 3D has a built-in counter for clock time. The 
wave patterns generated by its simulator will be updated at every 0.1 second interval of clock time. Unity 3D 
will start rendering its image as soon as it receives the buoy profile data. The process will generate a 
synchronized and smooth motion of waves and the buoy at 10 frames per second in the virtual environment. 

Lab VIEW® is the software responsible for controlling the hardware of the WEC mechanism via a C- 
Rio 9012 controller. It was programmed to accept the data from the MATLAB® simulation code and then use 
this data to actuate the hardware in real time. The controller was programmed to listen to a pre -identified port in 
the network at all times. In addition to the controller, a WEC mechanism was built, as shown in Fig. 2, in order 
to convert the mechanical energy into electricity and, thus, produce a hybrid WEC VR system with enhanced 
simulation capabilities. When the data sent by MATLAB® through UDP has arrived at the port, the C-Rio 
immediately seizes it. The buoy position data is stored in the C-Rio memory first, and then is utilized to 
calculate and transmit the right signal to the motor controller to replicate the heaving buoy motion provided by 
the simulation in MATLAB®. Since Lab VIEW® requires at least 20 data points of the buoy motion profile 
initially in order to calculate the contour that the stepper motor will follow, there is a slight time lag between the 
instant data is received and the instant the WEC mechanism starts to move. 

The WEC mechanism consists of a crane, a magnetic rod, a coil, an LED board, and a spring at the top. 
The crane, which is about 40 inches in height, was built to provide the mounting places and hold all the other 
components in place, as shown in Fig. 3. While the 6-inch-long magnetic rod travels up and down in sync with 
the buoy motion, the coil remains stationary. The 18-gauge diameter coil collects the change in flux created by 
the moving magnetic rod and converts it to current in order to light the LED on the circuit board. The spring at 
the top of the crane helps to hold the rod in place and allows it to move only in the desired direction. The motion 
of the magnetic rod is adjusted proportionally in the WEC mechanism so as to satisfy the limits of the power of 
the motor and the distance of the rail so that the heaving motion of the buoy is faithfully replicated in the 
hardware. For the selected wind speed and wave patterns selected in this study, no such adjustment is needed as 
the rail length can cover the vertical traveling distance of the buoy. 




Figure 2. Overview of hardware set-up for the WEC mechanism. 




Figure 3. Details of the crane used in the WEC mechanism. 
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V. DEMONSTRATION OF THE WAVE ENERGY CONVERSION PROCESS 

In this section, four cases with different wind speeds, 2, 4, 6 and 8 m/s, respectively, are simulated in 
order to demonstrate the working principle of the WEC process under various sea conditions. The wave profile 
is generated in the physics-based simulator using the unified directional spectrum from eq. (18) and assuming 
that the aforementioned wind speeds occur at 10 m above the sea surface and are sustained over a fetch of 100 
km; these conditions are typical of fully-developed seas. A half-submerged spherical buoy with radius of 1 m is 
utilized in the simulations. As the wind speed increases, the average wave amplitude increases as shown in Fig. 
4 (a) - (d). The corresponding motion of the buoy is also shown in these figures. The amplification of the 
motion of a WEC device is achieved when it operates at near-resonance conditions. A necessary condition for 
that is the buoy velocity be in phase with the wave excitation [25]. In regular waves, this can be achieved when 
the damping of the power-take-off (PTO) mechanism is equal to the radiation damping of the wave frequency 
[26]. When the WEC device is operating in irregular waves, the control of the buoy motion is typically 
performed by adjusting the damping of the power-take-off (PTO) mechanism either in a pseudo-continuous 
manner or via discrete control, e.g., latching [27]. In the current simulations, a simpler passive control approach 
is employed: the damping of the PTO mechanism is preset to match the radiation damping at the peak spectral 
frequency. Even though in this way the energy absorption in irregular waves is suboptimal, it still suffices to 
demonstrate the WEC process, without the requirement for prediction of the wave excitation [26]. 
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Figure 4. Buoy motion and wave profile (in meters) at different wind speeds. 

Next, the ocean waves are captured in the Unity 3D game scene, to show the working process of the 
wave energy conversion. The wind speeds in this virtual environment feature different ocean conditions from 
calm and glassy to rough seas as shown in Fig. 5 (a) - (d). As already mentioned, the WEC converts the heaving 
motion of the buoy to electric current. The time series of the electric current is collected using the WEC 
mechanism shown in Fig. 2 for each of the four cases. The corresponding power spectrum profiles collected 
from the WEC mechanism are plotted in Fig. 6 (a) - (d). It can be easily observed that there is no dominant 
frequency at which the electric current is generated. The root mean squared electrical power can be calculated 
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Apparently, while the difference in power between the 2 and 4 m/s wind speed cases is not significant, the two 
latter cases show an increase of 22.3% and 46.2%, respectively, compared to the 2 m/s wind speed case. The theoretical 
value of the root mean squared absorbed mechanical power of the buoy, i.e., neglecting friction and power take-off 
mechanism losses, obtained from eq. (28) and the corresponding electrical power for each of the four cases are presented in 
Fig. 7. This ratio is indicative of the WEC absorption efficiency. Even though the amount of absorbed power is greater at 
higher wind speeds, the efficiency drops mainly as a result of available hardware limitations, i.e., the limited useful range of 
change in magnetic flux of the utilized coil shown in Fig. 2. 




(c) v = 6 m/s (d) v = 8 m/s 

Figure 5. Ocean surface in Unity3D at different wind speeds. 
VI. CONCLUSIONS AND DIRECTIONS FOR FUTURE RESEARCH 

The development of a hybrid virtual reality (VR) software that simulates the wave energy conversion (WEC) 
process of a heaving buoy is described in this article. This VR system comprises three main function blocks: a physics-based 
simulator in MATLAB®, a game studio in Unity 3D, and a hardware manipulator in Lab VIEW®. The VR software is also 
combined with a WEC mechanism in order to demonstrate the conversion of mechanical energy into electricity. The 
software was tested and validated in irregular seas. The absorbed wave energy was calculated and the efficiency of the 
conversion process was quantified. The experiments demonstrated that the developed VR system can be utilized as both a 
research and an educational tool. Furthermore, the versatility of the employed wave spectrum regarding the modeling of 
irregular waves and the modularity of the VR system tools allow for easy implementation of a wide range of WEC 
simulation scenarios. 

A number of improvements that would enhance the effectiveness of the hybrid VR WEC simulation 
system are currently considered. First, as already mentioned, the utilization of a passive control strategy for the 
heaving buoy motion is far from optimal, thus, a more effective control method will be implemented in the near 
future. In addition to this, the utilization of an electromagnetic coil that would allow for increased conversion 
efficiency and easier manipulation of the amount of generated electricity is currently under way. Finally, the 
code developed for the VR system will be optimized in order to minimize the time lag in the communication 
between its three function blocks. 
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